Limited dependent variables

The problem

• Many important outcomes in the social sciences are non-continuous

  • Binary or categorical outcomes: vote choice, choosing to ride a bike instead of driving, or winning an election

  • Counts: number of trade agreements, number of gun deaths in a city, or number of doctor visits in a year

  • Durations: length of a military conflict, the longevity of a coalition, or the time spent on unemployment

The problem

Fitting a line to this kind of outcome often produces nonsensical results: predictions below the minimum or above the maximum and slopes that have impossible increments. Moreover, the functional form is likely incorrect even where the predictions are possible because the slope depends on the current expected value of the dependent variable.

The problem

And the resulting standard errors will be non-normal and heteroskedastic by definition:

NAVCO data

The NAVCO project collects data on “maximalist campaigns”, which are defined as sustained efforts to remove a government, attain territorial independence, or make fundamental changes to an existing constitutional order. A key question motivating the project is “what sorts of tactics are effective at creating major political change?”

navco|>select(LOCATION, CAMPAIGN, BYEAR, SUCCESS, VIOL, DURATION, REGAID, STATESUP, goal_type)|>head()

Some options other than logit/probit

Solutions: contingency table methods

• For models with a single categorical predictor and a limited number of outcome categories, contingency table methods can give you a valid test of statistical significance

  • The \(\chi^2\) and Fisher’s Exact Test will work for categorical by categorical relationships.

  • Spearman’s \(\rho\) , Kendall’s \(\tau\), and Somer’s \(D\) can test ordered relationships.

Solutions: linear probability model

Solutions: machine learning methods

Many machine learning methods can handle all sorts of weird data structures and many of them perform better than logit/probit models when predicting binary outcomes.

However! Machine learning means optimizing for prediction rather than hypothesis testing.

Logit and probit models

(The following sections will deal with models for binary outcomes, but the general intuition applies to other kinds of discrete choices)

The Latent Variable Intuition

The Logistic Function

Changing the “intercept” term here will shift the curve to the left or the right

The Logistic Function

Changing values of the slope will make the line steeper (especially around the middle)

Maximum likelihood

• OLS tries to minimize the sum of squared errors, but we can’t actually infer that here

• Instead, maximum likelihood methods attempt to identify the function that maximizes the likelihood of the observed data.

• There’s no “closed form solution”, so MLE methods use optimization algorithms to fit models

Maximum likelihood

\[ \ln L(\beta_0,\beta_1) = \sum_{i=1}^N \{ y_i \ln p(x_i; \beta_0,\beta_1) + (1-y_i) \ln [1-p(x_i;\beta_0,\beta_1)] \]

Maximum Likelihood

set.seed(400) 
N <-5000 
alpha <- 0
beta <- 1
d<-tibble(
  X = rnorm(N), # random IVs 
  y_star = alpha + X * beta,  # a latent variable
  Y = as.numeric((y_star + rlogis(N)) >0) # 1s and 0s (This is what we can observe!)
)

# the negative log likelihood function for the logistic model: 
neg_log_likelihood = function(beta, x, y){
  theta = beta * x
  p <- exp(theta) / (1 + exp(theta))
  val = -sum(y * log(p) + (1-y)*log(1-p)) 
  return(val)
}

Maximum Likelihood

Note that the negative log likelihood is minimized at the the correct value of \(\beta_1\):

The only thing you really need to know here, though, is that this is the thing we’re optimizing instead of minimizing the RSS. Maximizing this likelihood is non-deterministic, so it has to be approximated using an iterative process. This means it can fail to converge in some situations (usually small samples or lots of covariates)

Logit vs. Probit

The only real distinction here is that logistic models assume a logistic error, and probit models assume a Gaussian error. The interpretation of the coefficients will be different, but once they’re converted to predicted probabilities and marginal effects they should return very similar results.

set.seed(100)
N<-1000
X<-rnorm(N, 0, 1)
ystar<-X * 2 
# the probit link
y_probit<-(ystar + rnorm(N))>0
# the logit link 
y_logit <-(ystar + rlogis(N))>0
probit<-glm(y_probit ~ X, family=binomial(link = probit))
logit<-glm(y_logit ~ X, family=binomial(link = logit))

Interpretation

The way we solve this is ultimately less important than how it changes interpretations. There’s a few key differences here when compared to OLS

Side note: uncontrolled heterogeneity matters here

set.seed(100)
N<-1000
X<-rnorm(N, 0, 1)
Z<- rnorm(N, 0 ,1)
ystar<-X * 4  + Z * 2
y_logit <- (ystar + rlogis(N))>0
model1<-glm(y_logit ~ X, family=binomial(link = logit))
model2<-glm(y_logit ~ X +Z, family=binomial(link = logit))

Applying the model

navco_logit<-glm(SUCCESS ~  goal_type, data=navco, family=binomial(link='logit'))

navco_probit<-glm(SUCCESS ~  goal_type, data=navco, family=binomial(link='probit'))

modelsummary(list(
  "navco LPM" =  navco_lpm,
  "navco logit" =navco_logit,
  "navco probit"= navco_probit
                  
                  ), 
              coef_rename=TRUE,
             coef_omit = "Intercept",
              estimate  = "{estimate}",  
             statistic = c("conf.int"),
             conf_level = .95,        
 note = "95% CI in brackets",
 gof_omit = 'F|RMSE|R2$|AIC|Log.Lik.',
 output ='kableExtra'
             )

Interpretation: logit intercept

For logit, the coefficients are \(ln \frac{p}{1-p}\) of Y=1. So exponentiation can give you some sense for the probabilities if you’re used to thinking in those terms. After exponentiation, the intercept is the log odds when everything = 0.

# exponentiated coefficent
exp(coef(navco_logit))

                (Intercept)      goal_typeRegime Change 
                  0.6000000                   1.4885496 
         goal_typeSecession goal_typeSelf determination 
                  0.1937985                   1.8115942 

#  Pr(Y=1) at viol = 0 
pr<-mean(navco$SUCCESS[which(navco$VIOL==0)])
pr

[1] 0.5731707

# expressed as an odds ratio = the exponentiated intercept
orat<-pr/(1-pr)
orat

[1] 1.342857

Interpretation: logit coefficient

And the slope coefficient is the impact of a one unit change in X on the relative odds ratio compared to the baseline.

# exponentiated coefficent
exp(coef(navco_logit))

                (Intercept)      goal_typeRegime Change 
                  0.6000000                   1.4885496 
         goal_typeSecession goal_typeSelf determination 
                  0.1937985                   1.8115942 

#  Pr(Y=1) at viol = 1
pr_viol<-mean(navco$SUCCESS[which(navco$VIOL==1)])

pr_viol

[1] 0.2978723

# the odds ratio
orat_viol<-pr_viol/(1-pr_viol)

# the relative odds ratio compared to the baseline
orat_viol/orat

[1] 0.3159252

(Probit coefficients don’t work this way, which is probably why they’re less popular)

Interpretation: predicted probabilities and marginal effects

Rather than try to interpret logit/probit coefficients in their untransformed states, its almost always better to convert them to predicted probabilities and marginal effects at some meaningful levels of the coefficients.

Predictions

And automatically generate a nice plot for the predicted results:

plot(preds, 
     show_data = TRUE, jitter=.03)  # plot some relevant data points